#include <vector>
#include <iostream>
#include <algorithm>	//for find, min_element, for_each
#include <complex>
#include <iomanip>	//for setw
using namespace std;


//void find_pivot(vector<double> &C, vector< vector<double> > &A, vector<double> &b)

struct pivot{
	int row;
	int col;
	bool not_found=0; 	//if true will indicate a pivot could not be found
};

template<typename type>
struct CAb{
	vector<type> C;
	vector< vector<double> > A;	//pretty sure I can make this real instead of complex... look at it again though since I don't know what the top level functions give me	
	vector<double> b;
	type b0;
};

template<typename type>
void print_matrix(vector< vector<type> > &matrix)
{
	cout<<endl;
	for(int r=0; r<matrix.size();r++)
	{
		for(int c=0; c<matrix[r].size();c++)
		{
			cout.setf(std::ios::left, std::ios::adjustfield);
			cout<<setw(30)<<matrix[r][c];
		}
		cout<<endl;
	}
	cout<<endl;		
}

template<typename type>
void print_vector(vector<type> &vec)
{
	for_each(vec.begin(), vec.end(), [](type value){cout<<value<<"\t";});
	cout<<endl;
}

void print_tableau(CAb<complex<double>> &data)
{
	vector<vector<string> > matrix;
	ostringstream stream;
	vector<string> str_version;
	for_each(data.C.begin(),data.C.end(),[&str_version,&stream](complex<double> value){stream<<value;str_version.push_back(stream.str());stream.str("");});
	stream<<data.b0;	
	str_version.push_back(stream.str());
	matrix.push_back(str_version);
	stream.str("");
	str_version.clear();
	
	int j=0;
	for_each(data.A.begin(), data.A.end(),[&](vector<double> vec)	//adds Matrix A and its b's to the string matrix
	{
		for_each(vec.begin(),vec.end(),[&](double value){stream<<value;str_version.push_back(stream.str());stream.str("");});
		stream<<data.b[j++];
		str_version.push_back(stream.str());
		matrix.push_back(str_version);
		stream.str("");
		str_version.clear();
	});

	print_matrix(matrix);
	
}

bool check_real_less_than_zero(CAb<complex<double>> &data)	//checking C
{	
	//add tolerance code to round very small numbers to 0
	bool flag=0;
	for_each(data.C.begin(), data.C.end(), [&flag](complex<double> value){if(real(value)<0) flag=1;});	
	return flag;
}

bool check_imag_less_than_zero(CAb<complex<double>> &data)	//checking C
{
	//add tolerance code to round very small numbers to 0
	bool flag=0;
	for_each(data.C.begin(), data.C.end(), [&flag](complex<double> value){if(imag(value)<0) flag=1;});	
	return flag;
}	

struct real_less{
	bool operator()(complex<double> &val1, complex<double> &val2) const { return real(val1)<real(val2);}
};

struct imag_less{
	bool operator()(complex<double> &val1, complex<double> &val2) const { return imag(val1)<imag(val2);}
};


struct min_pos_comp{
	template<typename type>
	bool operator()(type &val1, type &val2) const
	{ 
		if(real(val1)<0 && real(val2)>=0)
			return 0;
		if(real(val1)>=0 && real(val2)<0)
			return 1;
		return real(val1)<real(val2);
	}
};

//vector<complex<double>>::iterator index_of_lowest_real(CAb &data) //can use auto but this is more informative
int index_of_lowest_real(CAb<complex<double>> &data)
{
	vector< complex<double> >::iterator index;
	index = min_element(data.C.begin(), data.C.end(),real_less());
	return index-data.C.begin();	//converts the iterator to int
}

int index_of_lowest_imag(CAb<complex<double>> &data)	
{
	vector< complex<double> >::iterator index;
	index = min_element(data.C.begin(), data.C.end(),imag_less());
	return index-data.C.begin();	//converts the iterator to int
}

pivot find_pivot(CAb<complex<double>> &data) 	//should return a triple (row,col,isnan)
{
	pivot pivot_point;
	if(check_real_less_than_zero(data)==1 || check_imag_less_than_zero(data)==1)	//check if this statement is necessary
	{
		int col;
		if(check_imag_less_than_zero(data)==1)	//imaginary first since it represents the coefficient of M and is much larger
		{
			col=index_of_lowest_imag(data);
		}
		else
		{
			col=index_of_lowest_real(data);	//case is ensured if the first doesn't occur. the top level if statement ensures this
			if (imag(data.C[col])!=0)		//if it is -a+Mb it doesn't attempt to pivot around this point. This is not a potential pivot since -a+Mb is a positive number
			{
				pivot_point.not_found=1;
			}
		}
		pivot_point.col=col;
		vector<double> temp_vec(data.A.size());
		for(int h=0; h<data.A.size();h++)
		{
			temp_vec[h]=real(data.b[h])/real(data.A[h][col]);	//A is in standard form so all entries are real... only possible none real entry is b[0] which we don't need			
		}	
		
		int row=min_element(temp_vec.begin(),temp_vec.end(),min_pos_comp())-temp_vec.begin();
		if (data.A[row][col]<=0)	//logic error changed. see simplex.cpp for former
			pivot_point.not_found=1;
		else
			pivot_point.row=row;		
		
	}
	else
	{
		pivot_point.not_found=1;
	}
	return pivot_point;
}


void basis_exchange(pivot piv,CAb<complex<double>> &data)		//column reduction
{
	int row=piv.row;
	int col=piv.col;
	cout<<"pivot:"<<"index:["<<row<<"]["<<col<<"] value: "<<data.A[row][col]<< endl;
	double denom=data.A[row][col]; //if i don't use this then for_each doesn't work correctly... strange error where data.A[row][col] value changes ... should be the same 
	for_each(data.A[row].begin(),data.A[row].end(),[&denom](double &value){value/=denom;});
	data.b[row]/=denom; 	//row+1 because i added an element b to represent the right side of the optimization function
				//consider instead adding an extra variable to store the right side of the objective function... that way b can be a vector of doubles instead of complex<double>'s
	for(int i=0;i<data.A.size();i++)
	{
		if(i!=row)
		{
			double temp=(-1)*data.A[i][col];	//this has to be here otherwise some columns don't change at all
			for(int j=0;j<data.A[i].size();j++)
			{
				data.A[i][j]=temp*data.A[row][j]+data.A[i][j];
				//newA(i,:)=((-1)*A(i,col)).*newA(row,:) + A(i,:);
				
			}
			//newb(i+1)=((-1)*A(i,col)).*newb(row+1) + b(i+1);
			data.b[i]=temp*data.b[row]+data.b[i];
		}
	}
	double neg=-1;			//this is necessary, above something similar occurs but the compiler converts (int)*(double) into double ... (int)*(complex<double>)!= complex<double>
	complex<double> temp2=neg*data.C[col];		//this is also necessary
	cout<<data.C[col]<<" vs "<<temp2<<endl;
	for(int j=0;j<data.C.size();j++)
	{
		//cout<<" temp2: "<<temp2<<"data.A[row][j]: "<<data.A[row][j]<<" data.C[j]: "<<data.C[j]<<endl;
		data.C[j]=temp2*data.A[row][j]+data.C[j];
		
	}
	data.b0=temp2*data.b[row]+data.b0;
}	


bool all_positive_coefficients(CAb<complex<double>> &data)
{
	int flag=1; //true that all coefficients are positive
	//for_each(data.C.begin(),data.C.end(),[&flag](complex<double> value) 	//can't use break within for_each
	for(int k=0;k<data.C.size();k++)	
	{
		//cout<<"imag(data.C[k]L "<<imag(data.C[k])<<endl;
		//cout<<"real(data.C[k])<0 "<<(real(data.C[k])<0)<<endl;
		if(imag(data.C[k])<0) 
		{	
			flag=0;	
			break;	//break loop
		} //imaginary part is just a representation not an actual imaginary number . it represents infinity 
		  // this is called big M method
		else if(imag(data.C[k])==0 && real(data.C[k])<0) //if imag part is positive then the number is positive 
		{
			flag=0;
			break;
		}
	}
	return flag;	
}

void solve(CAb< complex<double> > &data)
{
	print_tableau(data);
	pivot the_pivot=find_pivot(data);
	while(!the_pivot.not_found)
	{
		basis_exchange(the_pivot,data);
		the_pivot=find_pivot(data);
		print_tableau(data);
	}
	//add code here to indicate whether solution is optimal, infeasible, or unbounded
	//unbounded if the RHS of objective function has an imaginary number
	//infeasible if coefficients of objective function are not all positive
	//optimal if the coefficients of objective function are all positive and the RHS of objective function is real	

	bool check_positive=all_positive_coefficients(data);
	//cout<<"check_pos "<<check_positive<<endl;
	if(check_positive==1 && imag(data.b0)!=0)
	{
		cout<<"There is no feasible solution"<<endl;
	}
	else if(check_positive==1 && imag(data.b0)==0)
	{
		cout<<"An optimal solution exists"<<endl;
	}
	else if(check_positive==0)
	{
		cout<<"Problem is unbounded"<<endl;
	}
}


int main()
{
	//	int rows=3;
//	int cols=5;	
//	vector<vector<double> > matrix(rows, vector<double>(cols));
//	vector<vector<complex <double>>> A;	
	vector<vector<double>> A;	
	vector<double> r1{5,3,1,1,0,0,0};
	vector<double> r2{-5,6,15,0,1,0,0};
	vector<double> r3{2,1,1,0,0,-1,1};
	A.push_back(r1);
	A.push_back(r2);
	A.push_back(r3);

	vector<complex <double> > C{{-50,-2},{-75,-1},{-60,-1},0,0,{0,1},0};
	vector<double > b{144,240,80};	//the zero is added to the original vector b size of b=1+num_rows_of(A)... 1 is for C
	CAb<complex<double>> linprob;
	linprob.A=A;
	linprob.C=C;
	linprob.b=b; 
	linprob.b0={0,-80};	
	solve(linprob);	//problem is infeasible

	vector<vector<double>> A1;	
	vector<double> r11{3,-3,2,8,1,0,0};
	vector<double> r21{5,6,-4,-4,0,1,0};
	vector<double> r31{4,-2,1,3,0,0,1};
	A1.push_back(r11);
	A1.push_back(r21);
	A1.push_back(r31);

	vector<complex <double> > C1{-3,-1,-5,-4,0,0,0};
	vector<double> b1{50,40,20};	//the zero is added to the original vector b size of b=1+num_rows_of(A)... 1 is for C
	CAb<complex<double>> linprob1;
	linprob1.A=A1;
	linprob1.C=C1;
	linprob1.b=b1; 
	linprob1.b0=0;
	//solve(linprob1);	//problem is unbounded
	
	vector<vector<double>> A2;	
	vector<double> q1{1,2,1,0,0};
	vector<double> q2{3,1,0,1,0};
	vector<double> q3{0,1,0,0,1};
	A2.push_back(q1);
	A2.push_back(q2);
	A2.push_back(q3);

	vector<complex <double> > C2{-2,-3,0,0,0};
	vector<double> b2{10,15,4};	//the zero is added to the original vector b size of b=1+num_rows_of(A)... 1 is for C
	CAb<complex<double>> linprob2;
	linprob2.A=A2;
	linprob2.C=C2;
	linprob2.b=b2; 
	linprob2.b0=0;
	//solve(linprob2);	//has optimal solution
	
	//auto g=std::find(C.begin(),C.end(),<0);
	//print_matrix(linprob.A);
	//cout<<check_real_less_than_zero(linprob)<<endl;
	//cout<<check_imag_less_than_zero(linprob)<<endl;
	//print_vector(linprob.C);
	//cout<<index_of_lowest_real(linprob)<<endl;
	//cout<<index_of_lowest_imag(linprob)<<endl;

	//vector<double> tester{-1,-2,-3,-6,-100};
	//int index=min_element(tester.begin(),tester.end(),min_pos_comp())-tester.begin();
	//pivot check=find_pivot(linprob);
	//cout<<"pivot: "<<check.not_found<<":"<<check.row<<","<<check.col<<endl;
	//cout<<index_of_lowest_imag(linprob)<<endl;

//	for_each(linprob.A[1].begin(),linprob.A[1].end(),[](complex<double> &value){cout<<value<<" ";});
	//for_each(linprob.A[1].begin(),linprob.A[1].end(),[](complex<double> &value){value=value+1;});
	//print_tableau(linprob);
	//complex<double>g=1;
	//g=g/3;
	//pivot y;
	//y.row=2;
	//y.col=4;
	//basis_exchange(check,linprob);
	
	//cout<<linprob.A[1][2]<<endl;
	//for_each(linprob.A[1].begin(),linprob.A[1].end(),[&linprob](double &value){cout<<(linprob.A[1][2])<<endl;value/=linprob.A[1][2];});
	//for_each(linprob.A[1].begin(),linprob.A[1].end(),[&linprob](double &value){cout<<linprob.A[1][2]<<endl;value=1;});
	//print_matrix(linprob.A);
	//print_tableau(linprob);
}
